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3 . ABSTRACT 
< 

Radio waves propagating from distant pulsars in the interstellar medium 
(ISM) , are refracted by electron density inhomogeneities, so that the intensity of 
O observed pulses fluctuates with time. The observed pulse shapes are used to di- 
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agnose electron-density distribution in the ISM. The theory relating the observed 
pulse time-shapes to the electron-density correlation function has developed for 
30 years, however, two puzzles have remained. First, observational scaling of 
pulse broadening with the pulsar distance is anomalously strong; it is consistent 
with the standard model only when non-uniform statistics of electron fluctuations 
along the line of sight are assumed. Second, the observed pulse shapes are consis- 
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tent with the standard model only when the scattering material is concentrated 

52 ■ in a narrow slab between the pulsar and the Earth. 

C3 ■ 

We propose that both paradoxes are resolved at once if one assumes stationary 
and uniform, but non- Gaussian statistics of the electron-density distribution in 
the interstellar medium. Such statistics must be of Levy type, and the propagat- 
ing ray should exhibit a Levy flight rather than the Gaussian random walk implied 
by the standard model. We propose that a natural realization of such statistics 
may be provided by the interstellar medium with random electron-density dis- 
continuities. A Levy distribution has a divergent second moment, therefore, the 
standard approach based on the electron-density correlation function does not 
apply. We develop a theory of wave propagation in such a non-Gaussian random 
medium, and demonstrate its good agreement with observations. The qualitative 
introduction of the approach and the resolution of the anomalous-scaling paradox 
was presented earlier in [PRL 91, 131101 (2003); ApJ 584, 791 (2003)]. 



Subject headings: ISM: general — pulsars: general — scattering — MHD — turbulence 
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1. Introduction. 

Observations of pulsar signals provide a valuable tool for investigating electron density 
distribution in the interstellar medium. The pulsar intrinsic signal is narrow in time, being 
about few percent of the pulsar period. The observed signal (averaged over many periods 
of pulsation) is broad and asymmetric with a sharp rise and a slow decay. The pulse broad- 
ening is attributed to the random refraction the waves experience while propagating in the 
interstellar medium. 



1.1. Overview. 

The shapes of the signals, and their scalings with wavelength A and with the pulsar 
distance have been investigated observationally and analytically for 30 years. The standard 
theory of interstellar scintillations (described in detail below) assumes that the propagat- 
ing wave is refracted by random Gaussian gradients, thought to be a good approximation 
due to the central limit theorem (Tatarskii 1961). The symmetric Gaussian distribution is 
fully characterized by its second moment, and the standard theory of scintillations aimed 
at reconstructing this moment from observations. The contradictions of this theory with 
observations were noted in the early 1970's, although not many pulsars were investigated at 
that time to make a definitive conclusion. 

The signal shape is characterized by its time- width, r, estimated at the 1/e level. In 
1971, Sutton analyzed scalings of pulse broadening times, r, with the radio wavelength, A, 
and with dispersion measures, DM, corresponding to the pulsar's distance along the line 
of sight, DM = J dzn(z). Observational scaling for large dispersion measures, DM > 
30pc cm~ 3 , is close to r oc A 4 DM 4 , while the theory gave r oc \ 4 DM 2 . The recent ob- 
servational results are shown in Fig. (1). To overcome the difficulty with the anomalous 
DM-scaling, Sutton suggested that the interstellar turbulence was not statistically uniform 
along the lines of sight, so that the lines of sight for more distant pulsars intersected regions 
with stronger turbulence. 

The second paradox was encountered by Williamson (1972, 1973, 1974), who compared 
the observed shapes of the pulses with the shapes predicted by the standard theory. He 
obtained a surprising result that the model of continuous turbulent medium was not consis- 
tent with observations. Rather, the best agreement was given by the models where all the 
scattering material was concentrated either in a thin screen or in a slab covering only 1/4 
of the line of sight between the pulsar and the Earth, with nearly equal quality of fit to the 
data [see the discussion and Fig. 4 in (Williamson 1974)]. 
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These two assumptions that are necessary to reconcile the standard theory with the ob- 
servations may have physical grounds since the interstellar medium is not uniform and denser 
regions (with stronger turbulence) are encountered closer to the galactic center. However, in 
this paper we discuss a simple and perhaps more plausible physical explanation for the Sut- 
ton and Williamson paradoxes. According to this explanation, the Sutton and Williamson 
paradoxes are not the consequence of the non-uniform, large-scale galactic distribution of 
the electron density responsible for scattering. Rather, they reflect the universal properties 
of the microscopic structures of the density fluctuations. 



1.2. Levy flights. Mathematical background. 

In our recent work (Boldyrev & Gwinn 2003a,b, 2004) we proposed that the Sutton 
paradox is the evidence of statistically uniform but non- Gaussian electron distribution in 
the interstellar medium. We noted that the time-broadened pulses belong to those observa- 
tional objects that depend not on the moments of the electron-density distribution function, 
but on its full shape. The standard Gaussian approach did not recognize this fact because 
in the Gaussian case, knowledge of the second moment is equivalent to knowing the distri- 
bution function itself. However, for a general non-Gaussian electron-density distribution the 
second moment may formally diverge, and the electron-density correlation function may not 
characterize the pulse shape. 

Physically, this divergence means that the second moment of ray-angle deviation is 
dominated not by the bulk of its distribution function, but by the extremely far cut-offs at 
the tails of this function [as we will see in the example of Section 4]. In contrast, the observed 
signal shape is determined by the bulk of the distribution function and is not sensitive to 
its far-tail cut-offs. In other words, the electron-density correlation function and the shape 
of the observed pulsar signal provide different, complementary descriptions of the electron 
distribution in the ISM. 

We also note that in the Levy model, the pulse shapes may provide a unique opportunity 
to diagnose plasma fluctuations at very small scales, > 10 8 cm. These scales can be smaller 
than the Coulomb mean-free path in HII regions, and can be close to the ion gyro radius. 
Observations of pulse shapes may thus provide diagnostics for magnetized plasma fluctua- 
tions, which may help to elucidate the role of magnetic fields in interstellar turbulence, see 
also (Armstrong, Rickett, & Spangler 1995; Goldreich & Sridhar 1997). 

In the theory of scintillation, the quantity of interest is the transverse electron-density 
difference, averaged along the line of sight, An± = J dz[n(x. 1: z) — n(x 2 , z)] . A wave 
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propagating in the interstellar medium is scattered (refracted) by density inhomogeneities, 
and the statistics of the refraction angle is related to the statistics of density differences 
across the line of sight. 

In the approach proposed in (Boldyrev & Gwinn 2003a,b) we assumed that the distribu- 
tion function of An± was stationary and uniform, but non-Gaussian, and it had a power-law 
decay as |An_i_| — > oo. If this distribution does not have a second moment, then the sum of 
many density differences along the line of sight does not have the Gaussian distribution, i.e. 
the Central Limit Theorem does not hold. Instead, the limiting distribution, if it exists, must 
be the so-called Levy distribution, which is common in various random systems (Klafter, Zu- 
mofen & Shlesinger 1995). By construction, such a distribution is stable under convolution, 
i.e. the appropriately rescaled sum of two independent variables drawn from the same Levy 
distribution, has again the Levy distribution. The Gaussian distribution is a particular case 
of this more general stable distribution. 

The Fourier transform (the characteristic function) of an isotropic Levy probability 
density function, Pp(o), has the simple form, 

oo 

F(n)= J daPpia) exp(i(j, ■ a) = exp(-|/i| /3 ), (1) 

— oo 

where the Levy index (3 is a free parameter, < (3 < 2. This form is valid for an arbitrary 
dimensionality of the fluctuating vector a. As we explain below, in our case this vector 
is two-dimensional, a oc y J dz [n(z,x.i) — n(z, x 2 )], where x is the coordinate in the plane 
perpendicular to the line of sight, z, and y is the unit vector in the direction connecting the 
two points in this plane, y = y/|y|, where y = xi — x 2 . 

Equation (1) can be taken as the definition of the isotropic Levy distribution, and in the 
present paper we will consider only isotropic distributions. However, it is important to note 
that a stable distribution need not be isotropic [see, e.g., (Nolan 2002)]. In the simplest case, 
the anisotropic distribution is a copy of the isotropic distribution, rescaled along one axis; 
but far less symmetric Levy distributions are possible. Non-isotropic stable distributions 
might arise in the case of the magnetized interstellar medium, when wave scattering has one 
or several preferred directions. For example, anisotropic scattering may be characteristic 
of interstellar MHD turbulence, as described by Goldreich & Sridhar (1997); Lithwick & 
Goldreich (2001); Cho, Lazarian, & Vishniac (2000); Chandran & Backer (2002); Scalo & 
Elmegreen (2004), where the preferred direction is given by the large-scale magnetic field. 

The sum of iV Levy distributed variables scales as IX^ "! ~ N 1 ^, which becomes 
diffusion in the Gaussian limit (3 = 2. For f3 < 2, the probability distribution function has 
algebraic tails, Pp(a)a da ~ |o-|~ 1 ~ /3 da, for large and its moments, (|cr| < >), of the order 
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C > (3 diverge. In (Boldyrev & Gwinn 2003a,b), we considered smooth density fluctuations, 
a oc |xi — x 2 |, and the proposed Levy statistics of ray-angle deviation lead to the pulse- 
broadening scaling r oc \ 4 DM^ 2+ ^^ . This agrees with observations when (3 ~ 2/3, and 
provides a natural resolution to the Sutton paradox. 

1.3. The results of the paper. 

In this paper we consider scintillations produced by a general, non-smooth density fluc- 
tuations, a oc |xi — X2| a / 2 ; the limit of smooth density fluctuations corresponds to a = 2. 
More precisely, exponent a denotes the scaling of the density-difference PDF defined in for- 
mula (1), i.e. this PDF should have the form Pp(cr, y) ~ Pp(a /y a l 2 )/y a l 2 . In the non-smooth 
case, the effective Levy index that enters the expressions for angle deviations is different from 
j3, and is given by 7 = a/3/2, as is explained in formulae (9) and (17) below. 

As we noted above, in the Levy case the electron-density correlation function (or second 
order structure function) is contributed to by the far-tail cut-off of the electron-difference 
distribution function. We will denote the correlation-function scaling as (a 2 ) oc \y\ s , where 
the scaling exponent, S, is in general different and independent of the Levy distribution 
scaling exponents a and (3. The exponent 5 is related to the Fourier spectrum of electron 
density fluctuations, {\n(k)\ 2 )k 2 oc k~~ 5 . 

In the present paper we demonstrate that the proposed Levy statistics of ray-angle 
deviations provide a resolution for the Williamson paradox as well. For this purpose we, 
first, develop a theory of wave propagation in a Levy random medium by approximating 
the medium by a uniform series of scattering screens. Then, we provide a general method 
for constructing the pulse-broadening function for an arbitrary Levy index 7. And finally, 
we compare our results with the observational signals of large-dispersion-measure pulsars, 
recently published by Bhat et al (2004) and by Ramachandran, et al (1997). We obtain that 
the observational shapes agree well with the predictions of our theory for a uniform medium 
with the Levy index 7 ~ 2/3 to 1, while they are inconsistent with the Gaussian theory 
corresponding to the Levy index 7 = 2. 

We also discuss the effect of "over-scattering" which is inherent for the Levy scintilla- 
tions, because the time-signals have long algebraic tails that do not decay to zero during 
pulse periods. This effect is crucial for comparing analytical and observational data, as we 
show in Section 5. 

For illustrative purposes, we also present a particular model of the density distribution 
that produces the Levy statistics of scintillations. The proposed density distribution is 
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strongly spatially intermittent, and it can be visualized as follows. Imagine that the electron 
density is concentrated in separated regions that have sharp irregular boundaries. We make 
the simplest assumption that these boundaries are random, similar to randomly folded two- 
dimensional sheets or shocks. If one further assumes that the line of sight intersects the 
boundary at a random angle whose distribution is uniform in all directions, such a picture 
corresponds to the particular case of a Levy distribution of An, with a — 2, (5 — 1, the 
so-called Cauchy distribution. This distribution is rather distinct from the Gaussian one, 
and is close to the distribution that has been predicted in our theory. 



1.4. Related models. 

The importance of randomly oriented discontinuous objects that can be encountered 
across the line of sight has been emphasized in earlier theoretical work on scintillation, 
see e.g., (Lambert & Rickett 2000) and (Cordes & Lazio 2001). As noted by Lambert 
& Rickett (2000), sharp density discontinuities may characterize stellar wind boundaries, 
supernova shock fronts, boundaries of HII regions at the Stromgren radius, etc. They may 
also arise from supersonic turbulent motion. Rickett, Lyne, & Gupta (1997) proposed that 
similar highly intermittent density structures may be responsible for observed low-frequency 
modulations (fringes) of dynamic pulsar spectra. 

The approach of Lambert & Rickett (2000) utilized the non-Kolmogorov spectra of 
discontinuous density fluctuations. Cordes & Lazio (2001) considered wave scattering by 
confined or heavily modulated screens (such as disks, filaments, etc.), when the statistics 
of angular deviations were generally assumed Gaussian, with parameters varying along the 
screens. Both papers discussed important aspects of non-Kolmogorov and spatially inter- 
mittent electron-density distribution in the ISM, however, they implied the existence of the 
second moments of the scattering-angle distributions. Therefore, in the earlier considera- 
tions, the effects that we discuss in the present paper could not be discovered. 

As we demonstrate in section 4, sharp density discontinuities can, in fact, produce a 
non-Gaussian distribution of the scattering angle, whose second moment diverges. In our 
approach we exploit such intermittent density statistics in their full generality. To elucidate 
the universal role of Levy distribution, we keep our consideration as simple as possible, 
assuming statistically uniform and isotropic scattering screens (although generalizations for 
the non-isotropic and non-uniform cases are possible). The presence of strongly-scattering 
structures is naturally represented in our model by slowly decaying, power-law tails of the 
scattering-angle distributions. Our model has a simple physical interpretation and provides 
a practical way of calculating pulse shapes and pulse scalings. Most importantly, it naturally 
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resolves both the Sutton and the Williamson observational paradoxes, in a manner that is 
simpler and complementary to the standard Gaussian picture. 

The Levy-flight model may also be relevant to observations of extreme scattering events, 
such as those reported by Fiedler et al (1987); Wolszczan & Cordes (1987); Stinebring et 
al (2001). Such events have been investigated theoretically, as by Rickett, Lyne, & Gupta 
(1997); Romani, Blandford, & Cordes (1987). We note that these theoretical studies invoke 
uncommon incidents of scattering much larger than typical values to explain these events. In 
this sense, these events are consistent with the Levy model we develop here, which includes 
rare, large events in a statistically-stationary way. However, we do not yet know whether 
these phenomena might fit into a single model for interstellar scattering along with the 
common phenomena we discuss in this work. The answer might possibly depend on the 
details of the picture for turbulence in which the Levy model is realized. Accordingly, we do 
not discuss the extreme scattering events in the present work. 



To address the puzzles mentioned above we need to review the standard theory of 
interstellar scintillations. First, we note that in the interstellar plasma with typical electron 
density n ss 0.03 cm -3 , the electron plasma frequency is oo pe = (ATcne 2 /m e ) 1 ^ 2 10 4 s -1 . 
This frequency is much smaller than the typical observational frequency of 10 8 — 10 9 Hz, 
and, therefore, the propagating wave scatters only by a small angle on the scale of density 
inhomogeneities. The Fourier amplitude of electric field, E(cj,r) = J E(t, r) exp(— iuijdu, 
obeys the following equation: 



see, e.g., (Tatarskii 1961; Lee & Jokipii 1975a). We are not interested in polarization 
effects that are small in the considered approximation by a factor Q e /u ~ 10~ 8 , where 
Q e = eB/(m e c) is the electron-cyclotron frequency, therefore, we consider the scalar wave 
amplitude, E(cu,r). 

Equation (2) can be reduced further, using the so-called parabolic approximation (Tatarskii 
1961). Assuming that the wave propagates in the line-of-sight direction, z, we can sep- 
arate the quickly changing phase of the wave from the slowly changing wave amplitude, 
E(u,r) = exp(izu/c)$ u (z, x), where x is a coordinate perpendicular to z. The equation for 
the wave amplitude $ w reads: 



2. Directed-wave propagation in a random medium. 




(2) 
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where Aj_ is a two-dimensional Laplacian in the x plane, and r = e 2 /(m e c 2 ) is the classical 
radius of the electron. 

Following the approach of Uscinski (1974); Williamson (1975); Lee & Jokipii (1975a,b,c) 
we introduce the two-point function J(r 1 ,r 2 ,t) = $(ri, t)$*(r 2 , t), whose Fourier trans- 
form with respect to time is /^(ri,r 2 ) = l/y/2n J dLu^ UJ+ Q/ 2 (ri)^* uj _ n ^ 2 (r2). For coin- 
ciding coordinates, ri = r 2 , this function is the intensity of the radiation whose vari- 
ation with time we seek. To find this function we need first to solve the equation for 

-fV2( r 2)> wn i cn can be derived from Eq. (3). Assuming that ^ < w, 
we obtain the following equation, 

,8V 1 d 2 V 1 d 2 V 27rr 

% ~dz~ ~ ~2k ± 9xY + 2k 2 9xf + ~T~ V > { ) 

where we denote k = cu/c, An = n(^,xi) — n(z, x 2 ), and Ak = k\ — k 2 = Q/c. 

Equation (4) is hard to solve since n(z, x) is an unknown random function. The standard 
theory uses the fact that it takes many refraction events to appreciably deviate the ray 
trajectory, as we discussed in the introduction. Therefore, due to the central limit theorem, 
the ray deviation angle exhibits a Gaussian random walk. One can therefore assume that 
the density fluctuations are Gaussian with the specified second-order correlator ([n(zi,xi) — 
n(z 2 ,x 2 )] 2 ) = 2«(x 1 — x. 2 )S(zi — z 2 ), where short correlation length in z direction is the 
mathematical expression of the fact that the ray becomes appreciably deviated only when 
it travels the distance much larger than the the density correlation length (Tatarskii 1961; 
Lee & Jokipii 1975a,b,c; Rickett 1977, 1990). The correlation function n(y) is the two- 
dimensional Fourier transform of the spectrum of the density fluctuations. More precisely, if 
the three-dimensional spectrum of the density fluctuations is D{k) = (\nk\ 2 ), then n(y) is the 
Fourier transform of D(k±, k z = 0) with respect to kj_. When the density has a power-law 
spectrum, D(k)k 2 oc k~ 5 , with 1 < 5 < 2, then n(y) oc y 5 . 

As we mentioned in the introduction, the statistically uniform Gaussian model fails 
to reproduce the observational scaling of the the broadening time, r ~ A 4 DM 4 . Numer- 
ous attempts to reproduce this scaling by using different spectra of density fluctuations 
within this model have not been successful either. The purely shock-dominated density 
distribution has the spectrum (\nk\ 2 )k 2 oc k~ 2 , while the Kolmogorov turbulence has the 
spectrum (\rik\ 2 )k 2 oc /c™ 5 / 3 ; the difference in the spectral exponents is rather small to have a 
considerable consequence. For various important aspects of the standard Gaussian scintilla- 
tion theory we refer the reader to (Goodman & Narayan 1985; Blandford & Narayan 1985; 
Gwinn et al 1998; Lambert & Rickett 2000; Lithwick & Goldreich 2001). In other contexts, 
the theory of wave propagation in Gaussian random media was developed in Saul, Kardar 
& Read (1992); Jayannavar & Kumar (1982). 
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In the next section we formulate the problem by assuming that the refraction occurs in a 
series of discrete structures, refracting 'screens', that are statistically identical, independent, 
and placed uniformly along the line of sight. This setting is physically appealing since, as we 
have mentioned in the introduction, the refraction in the interstellar medium is consistent 
with the presence of spatially intermittent scattering structures, see e.g., (Lambert & Rickett 
2000; Cordes & Lazio 2001). Moreover, this is the simplest setting when the problem admits 
an exact analytic solution that allows us to treat both Gaussian and Levy cases on the same 
footing. 



3. Multi-screen model of scintillations. 

Let us assume that wave scattering occurs in a series of thin flat screens uniformly 
placed along the line of sight. Also, for simplicity assume that the incident wave is planar, 
although an analogous consideration can be made for spherical geometry as well. Each 
screen gives a contribution to the signal phase; if we denote the wave function just before 
the screen by V m (xi,x 2 ), then the wave function right after the screen will be given by 
\/ ou '(xi, x 2 ) = S , (xi,x 2 )y m (xi,x 2 ). The phase function, S, can be found from Eq. (4), 

S , (xi,x 2 ) = exp ^ 2 ^ rr ° J dz [n(z, x x ) - n(z, x 2 )] j , (5) 

where the integration is done over the thickness of the phase screen, /. Between the phase 
screens, the wave propagation is free. Effecting the change of variables, y = x x — x 2 , and 
x = x x + x 2 , we rewrite the Eq. (4) in this region, 

.W_Ak((fV_ d 2 V\ 1 d 2 V 
dz k 2 \ dy 2 <9x 2 / k <9x<9y ' 

In what follows we will be interested in the wave function, V, averaged over different 
realizations of the electron density in the phase screens. Due to space homogeneity, the 
averaged transfer function, S = (S'(xi,x 2 )), should depend only on the coordinate difference 
y. Since refraction affects only the y dependence of the wave function, we can assume 
that the averaged wave function is independent of x, (V(z, x 1; x 2 )) = U(z,y). Taking into 
account the x dependence would lead to more cumbersome formulas, although it would 
not qualitatively change the results. We therefore model the free propagation between the 
screens by 

.dU _Akd 2 U 
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To describe the scattering we need to specify the averaged phase function, S. In the 
Gaussian case, one can assume that the screen width, /, is of the order of the characteristic 
length of the density fluctuations, l , and therefore 1 , 

S = exp [-A 2 r V(y)] , (8) 

where A = 2n/k. Assuming that the main contribution to the scattering comes from the 
scales, y, much smaller than the density correlation length, we can expand n(y) K,o(y/lo) a , 
where k ~ Z (An ) 2 , and An is a typical amplitude (say, rms value) of density fluctua- 
tions. For example, if density fluctuations arise due to passive advection by the Kolmogorov 
turbulence, then a = 5/3. The smooth density field corresponds to a = 2. 

The situation is completely different in the Levy case, when the second-order moment 
of the density integral in (6) diverges. We however can assume that the density difference 
has some scaling form, An± = Jjj io [n(xi, z) — n(x 2 , z)\dz ~ l An (y/l ) a ^ 2 , and use the Levy 
formula (1) to write the most general expression for the averaged phase function S, 

S = exp (-|Aro/oAn | /3 | 2 ///o| Q/3/2 ) . (9) 

In the Gaussian limit, (3 = 2, this formula reduces to the previous result (8). Similar to the 
Gaussian case, the exponent a/2 has the meaning of the density-difference scaling with the 
point separation, however, it characterizes not the second moment (that does not exist), but 
the density-difference distribution function itself. 

More precisely, for a chosen point separation y, the projected density difference in the 
screen function (5), 

A<3> = Ar / dz [n(z, xi) - n(z, x 2 )] = Xr An±, (10) 
Jo 

should be drawn from a distribution function that depends only on the combination A$/|/ Q//2 . 
We also note that contrary to the Gaussian case, the quantity An entering equation (9), 
does not correspond to the rms value of density fluctuations, since (An 2 j_) does not exist in 
the Levy case. Rather, it denotes the typical width of the density-difference PDF, say the 
density difference at which this PDF decays by the factor 1/e compared to its maximum 
value. On the contrary, the rms value of An± would be given by far, non-Levy tails of this 
PDF, and would be much larger than An . 



1 The assumption I ~ Iq is not necessary. The only assumption required for our consideration is Iq < I, 
and the final results can be generalized for this case. This generalization would however require further 
assumptions about the scattering structures, which at the present level of understanding can hardly be 
justified. 
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Physically, the Levy screen function in the form (9) could correspond to the distribution 
of the phase integral, A<3>, when the integration path intersects a random density discontinu- 
ity, say a sharp boundary of an ionized region. Such a random boundary may naturally arise 
as a result of turbulent advection in an ionized region, and numerical simulations could test 
this intriguing possibility. Our suggestion is based on a result that a smooth, randomly ori- 
ented boundary reproduces the particular form of the Levy screen function (9) corresponding 
to P — 1, a — 2. We will derive this result in the next section. In general, particular values 
of a and (3 depend on the physical realization of scattering screens, which is not completely 
understood and cannot be specified in this work. For our present purposes, we simply assume 
the general form of the screen function (9) and provide a method for calculating the shape 
and the scaling of the received signal. 

Let us consider a series of such screens placed along the line of sight, z, perpendicular to 
it, such that the distance between two adjacent screens is zq. Number the screens from the 
source to the observer, and consider the mth screen. Denote as U™(y m ) the wave function 
just before the wave passes the screen, and L^ 1 * = SU™ this function just after the screen. 
Between the screens to— 1 and to, the propagation is free, therefore, we obtain from equations 
(4) and (7), 

We can write the wave function after the wave has passed N screens, by iterating this formula 
N times, and by using the expression for the screen function, (9), 

/ N N-l \ 

UftiVN) = [-inAAk] J exp [iAAk]' 1 £ (Ay m ) 2 - B°W £ \y m \^' 2 x 

\ m=l m=l / 

U (y )d 2 y N _ 1 ,...,d 2 y 1 d 2 y (12) 

where we introduced the short-hand notation: A = 2z /k 2 , B = \\r loAno\ 2 ^ a /l , and 
Ay m y m Ym-i- 

To do the integrals, we need first to make a simple transformation of formula (12). We 
will substitute the following identity, 

exp [(Ay m ) 2 /(iAAk)] = -irrAAk J exp (i£ m • Ay m - iAAkfa) d 2 U, (13) 

for each Ay m , into formula (12). Our ultimate goal is to find the time dependence of the 
pulse intensity at the Nth. screen, i.e., I(t) = J Un{un = 0, Ak) exp(iAk ct)dAk c. We 
therefore change the order of integration and do the A/c-integral first. As the result we get: 

/ N-l N \ 

I(t) = / exp -B^l 2 £ \y m \»V 2 - i £ £ m • Ay m U (y ) x 

\ m=l m=l / 
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8 (t ~ A £ &/c) rf 2 ^ • • • rf 2 6 <*W-i • • • d 2 y . (14) 

Introducing the non-dimensional variables y m = y m B and £ m = £ m /B, we rewrite formula 
(14) as follows: 

/ n-i n _ \ 

l(t) = Jexp - E |y m | a/3/2 -^EL-AyJ C/o(yo/S) x 

\ m=l m=l / 

s(t-*£ E x e) ^ • • • d% d 2 y N _, . . . d 2 y . (15) 

In the rest of this section we will use only the non-dimensional variables, and will omit the 
tilde signs. 

Now we are ready to do the y-part of the integral. For this we make a simple re- 
arrangement in the sum ^(y t - y ) + . . . + £ N (y N - yjv-i) = -jfofi + Vn£n + - 



£2) + • • • + Vn-i(£n-i — £n)- Recalling that y^ should be set to zero, we are left with 
Em=i£m ' A y m = -Vo£i ~ E™=J Ym • Af m+1 , where A£ m =£ m - £ m ~i- To complete the 
^-integration, we simply use the Levy formula (1), and get: 

r ( ab 2 N \ 

I(t) = J Si t - — Cn J Pam(^N) ■ ■ ■ P am (^ 2 )U (^B)d^ N . . . d%, (16) 



m=2 



where Uo(0 is the Fourier transform of Uo(y), P is the Levy distribution function introduced 
in (1). The last step is to note that d 2 ^^ ■ ■ ■ d 2 ^ — d 2 A^ . . . d 2 A^ 2 , which allows us to give 
a simple interpretation to formula (16). 

Before discussing this interpretation, we note the quantum mechanical analogy of our 
approach. If we could assume that B a/3 ^ 2 oc Zq, and zq oc 1/N in the limit when the 
number of screens increases, N — > 00, but the line-of-sight length is constant (which would 
be natural for the Gaussian case), the integral that we calculated in (12) would formally 
become the Feynman-Kac path integral for the solution of the Schrodinger equation with 
the potential oc i|y| Q/3//2 . The transition from (12) to (16) given by formula (13) is the 
transition from the coordinate to the momentum representation of this integral. 

Naturally, in the Gaussian case one can find the shape of I(t) either by doing the 
multiple integral in (16) or by solving the corresponding Schrodinger equation. The first 
approach was essentially adapted in (Uscinski 1974; Williamson 1975), the second one in 
(Lee & Jokipii 1975b). In our Levy model we cannot assume that B af3 ^ 2 ~ 1/N, therefore, 
we have to work with the general expressions (12)-(14). Although the general integral (12) 
is quite complicated, the quantity of interest, I(t), which is calculated with the aid of this 
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integral, has a quite simple meaning, and its probabilistic interpretation (16) can be easily 
understood. 

The 5-function in formula (16) means that the signal intensity I(t) is the probability 
density function of the time delay, r = AB 2 J2m=i £m/ c - The variable £ m is proportional 
to the deviation angle of the ray path from the z-axis, £ m oc m . [This is clear from (13) 
since £ m oc Ay m /zo, and this will also be the case in the geometric-optics analysis of Sec- 
tion 4.] This angle is the sum of elementary angle deviations caused by each phase screen, 
viz £ m = £i + Y^iL 2 The propagation time between two neighboring screens exceeds the 
propagation time along the z-axis by the amount At = (zq/c)[1 — cos# m ] oc and to get 
the total time delay, r, we need to sum up these individual time delays. 

Formula (16) teaches us that each angle increment, A£ s , has a Levy distribution with 
the index 7 = a/3/2, which provides a practical way of constructing the pulse broadening 
function I(t) without doing the multiple integrals in (16). This is a Monte-Carlo-type method 
for multiple integration. Indeed, to find the effect caused by the interstellar medium, we can 
neglect the intrinsic pulse shape, i.e., we can assume that the initial angle of the planar wave 
is zero, £ x = 0. Then the function I(t) is the probability density function of the variable 



where all the variables A£ s are distributed independently and identically according to the 
Levy law with the index 7 = a/3/2. The same expression was obtained in our previous work 
for the case a — 2, although without a detailed derivation (Boldyrev & Gwinn 2003a,b, 
2004). We remind that the scaling exponent a corresponds to the spatial scaling of the 
density field, as is defined in Eq. (9). The Gaussian case corresponding to a = 2, f3 = 2 
was originally considered by Williamson (1972). The Gaussian case corresponding to the 
Kolmogorov density fluctuations, a = 5/3, (3 = 2 was considered by Lee & Jokipii (1975c); 
however, the central part of the signal I(t) predicted by their model did not differ much from 
the Williamson case. 

We do not know the analytic expression for the probability density function of r. How- 
ever, it can be calculated numerically, by the method described in (Chambers, Mallows, & 
Stuck 1976; Nolan 2002). The generated pulse shapes for different numbers of screens, N, 
seem to be universal when appropriately rescaled. We will construct such shapes and will 
compare them with observations in Section 5. 

The scaling of the time delay r, i.e. the approximate scaling of its probability density 
function, can be found from the following consideration. Using the expression for the scaling 
of the sum of Levy distributed variables, | ^ N A£| ~ N 1 ^ (recall that the Central Limit 
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Theorem does not hold), we derive r oc A 2+4 / a A^^ 2+7 ^ 7 , where the number of screens is 
proportional to the distance to the pulsar, N oc DM. The observational scaling, r oc X 4 DM 4 , 
is reproduced for a ~ 2, and (3 ~ 2/3, which shows that ray-deviating density fluctuations 
do not have a Gaussian distribution but rather a Levy distribution with the index (3 ~ 2/3. 
This is the important result of the Levy model of scintillations, and it resolves the Sutton 
paradox. 

In the Section 5, we will demonstrate that the Williamson paradox is resolved natu- 
rally in the Levy picture of scintillations as well. But before that, in the next section we 
would like to present a physical model of scintillations that can be easily visualized and that 
demonstrates the essence of our approach. 



The intensity function I(t) can be calculated for the phase-screen function (9), once the 
parameters a and (3 have been specified from the physics of the problem. As an important 
example, let us consider a particular case of Levy phase screens, which we shall call Cauchy 
screens, when the values of these parameters can be derived directly. Let us assume that the 
interstellar medium is filled with separated ionized regions such that all of them have the 
same electron-density contrast with the ambient medium, An . We also assume that these 
regions are randomly shaped and have sharp boundaries. 

A propagating ray is then refracted by these random boundaries due to the Snell's law, 
771 sin(#i) = 772 sin(0 2 ), where 6 is the angle between the ray and the normal to the boundary, 
and i] = (1 — cjpg/cu 2 ) 1 / 2 m 1 — cj 2 e /(2u; 2 ) = 1 — (27r) _1 A 2 r n is the refraction index. As 
we explained in the introduction, the interstellar parameters ensure that ujpjuj 2 <C 1, and 
each refraction event results in a very small angle deviation, 56 = 61 — 62- Expanding 
sin(0i) = sin(# 2 + 56) through the first order in 56, we get from the Snell's law 



where An = n\ — n 2 . The distribution of tan(0) can be found assuming (somewhat artifi- 
cially) that the boundary normal is uniformly distributed over all directions in a half-space, 
at the point where it intersects the line of sight. Denote a = tan(#)n, where n is a unit vector 
indicating the direction of the refraction in the plane perpendicular to the ray propagation. 
Then the distribution of a is given in the polar coordinates by: 



4. Scintillations caused by random density discontinuities 



56 = (2n)- 1 \ 2 r An tan(6), 



(18) 



P(a) dad<j> 



1 



da d(f). 



(19) 



27T (1 + (7 2 ) 3 /2 
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Quite remarkably, this distribution is the two-dimensional Levy distribution with the index 
[3 = 1, also known as the Cauchy distribution. Its Fourier transform follows from formula (1), 
F(n) = exp(— \fj,\), we leave the derivation of this result to the reader as an instructive 
exercise. 

To apply this result we need to calculate the screen function (5). It can be done if we 
note that for the plane-like density discontinuity, the screen phase (10) is easily calculated, 
A$ = Ar / dz[n(z,yLi)—n(z,yL2)] = \r An (cr-y). Since a has the Cauchy distribution (19), 
the averaged screen function can be found with the aid of formula (1), 



and we recover formula (9) with / = l , and with a — 2, (3 — 1; the parameters are close to 
those that we have proposed for the interstellar scintillation on the observational grounds, 
i.e., a w 2, (3 xs 2/3. We expect that the observational scaling (3 ^ 2/3 would correspond to 
more complicated, realistic structure of the boundaries, possibly, with the fractal dimension 
larger than 2. Interstellar turbulence can indeed corrugate two-dimensional structures to 
make them have higher-than-two fractal dimensions, as is seen in numerical results and 
in observations of molecular clouds (Boldyrev 2002; Boldyrev, Nordlund & Padoan 2002; 
Kritsuk & Norman 2004; Elmegreen & Falgarone 1996); see also (Meneveau & Sreenivasan 
1990; Constantin, Procaccia, & Sreenivasan 1991; Kraichnan 1991). We expect this effect to 
have implications for interstellar scintillations, and we plan to consider it in the future work. 

An attentive reader has probably noticed that Eq. (18) for the angle deviation cannot be 
valid for the arbitrarily large tan(0). Indeed, the refraction angle 62 cannot exceed the critical 
angle, given by sin(# c ) = 771/772, which is obtained when 9\ = n/2, where we assume 772 > 771. 
To obtain the general formula for the angle deviation we need to expand sin^) = sin(92 + 59) 
in the Snell's law up to the second order in 59, which gives 



The distribution of 59 thus has a cut-off at 59 c = (A 2 r Ano/7r) 1 ^ 2 . For example, if the 
electron density fluctuations are An ~ 10~ 2 cm~ 3 , we have (27r) _1 A 2 roA77, ~ 10~ 12 1, 
which implies that formula (18) holds in the broad range of scales. Indeed, the distribution of 
a = 2ir59 I (\ 2 r An ) follows the Cauchy law (19), therefore the distribution function of 59 is 
strongly peaked at small 59. This function has a maximum at 59 ~ (27r) _1 A 2 roAn ~ 10~ 12 , 
and a long asymptotic power-law tail, P(59) oc l/(59) 2 , that spans about six orders of 
magnitude before it reaches the cut-off at 59 c ~ 10~ 6 . 

This consideration provides an illustration to the main idea of our approach. The second 
moment of the distribution of 59 depends not on the shape of the distribution function (19), 




(20) 




(21) 
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but on the cut-off value of its asymptotic power-law tail, 56 c . On the contrary, the shape 
of the pulsar signal is determined by the full shape of the distribution function (19), and is 
practically independent of the large cut-off value 59 c . 

5. Comparison with observations. 

In Fig. (2) we present the pulse-broadening function, I(t) = {5(t — r)) which is the 
distribution function of the time delay, r, given by formula (17), for the standard Gaussian 
model, a — 2, (5 — 2, for the Levy-Cauchy model, a — 2, (5 — 1, and for the Levy model with 
a — 2, /3 — 2/3. Note that, in comparison with the pulse through a Gaussian medium, the 
Levy models predict relatively greater intensity at short delays from the arrival time for the 
original pulse, as well as power-law tails at long delays. These tails result from rare, large 
occurrences of scattering. 

We generated the two-dimensional isotropic Levy distribution of A£ s using the methods 
of Feller (1971); Chambers, Mallows, & Stuck (1976); Nolan (2002). Let X = (X U X 2 ) be a 
Gaussian random vector, i.e. both components X\ and X 2 are independent, identically dis- 
tributed, one-dimensional Gaussian variables. Let Y be a completely skewed (i.e. positive), 
independent Levy distributed variable corresponding to the Levy index 7/2. This variable 
can be numerically generated using the method by Chambers, Mallows, & Stuck (1976), see 
also (Nolan 2002). Then the variable A£ = X\/F has an isotropic two-dimensional Levy 
distribution with the index 7; the proof can be found, e.g., in (Feller 1971). To generate an 
anisotropic Levy distribution with uniform index 7 (as may be relevant for the magnetized 
interstellar medium) one may choose an anisotropic Gaussian vector X, and then proceed 
as above. In our present work we consider only isotropic distributions. 

We compared the predicted pulse shapes with the data recently published in (Bhat 
et al 2004) and in (Ramachandran, et al 1997). Out of the 76 pulsars analyzed by Bhat 
et al (2004), we considered those with broad signals, in order to minimize the intrinsic- 
pulse effects. Good examples are provided by the four pulsars, P1849+0127 at 430 MHz 
(DM=214.4), J1852+0031 at 1175 MHz (DM=680.0), J1905+0709 at 430 MHz (DM=269.0), 
and P1916+0844 at 430 MHz (DM=339.0), where higher-frequency observations indicate 
that the intrinsic signal may be narrow. The dispersion measures are given in the units 
of pc ■ cm -3 . Shapes of all of the pulses have a characteristic sharp rise and a narrow, 
pointed apex, which is inconsistent with the Gaussian model, in exact agreement with the 
observations by Williamson (1974) [see Fig. (4) in Williamson's paper]. On the contrary, 
our Levy model provides a good fit to such shapes. For another similar comparison we used 
the pulse shape of J1848-0123 (DM=159.1) observed at 327 MHz by Ramachandran, et al 
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(1997); and we present this comparison in Figs. 4,5 as the illustrative example. 

When making the comparisons, we tried to match the central parts of the observational 
and analytical curves. However, this was possible only when we shifted the zero level of the 
analytical curves down, as is shown by the thin solid lines in Fig. (4). This may be natural 
for "over-scattered" profiles, i.e., when the signal has a long tail that does not approach 
zero during the pulse period. This is exactly the case for our Levy model. Besides, other 
observational effects, including noise from the sky and the telescope itself, usually prevent 
one from determining the baseline unambiguously. 

We note that for a precise comparison, we must convolve the calculated pulse-broadening 
function with the narrow intrinsic shape of the pulse, and with the response function of the 
receiving system Ramachandran, et al (1997); Bhat, Cordes & Chatterjee (2003). This will 
lead to slight broadening of the predicted pulse profile, making the Gaussian model even less 
consistent with the observations. 

The effects of over-scattering on baseline have not been considered observationally. This 
effect is important for comparison of our theory with observations, and we plan to investigate 
it in future work. Interferometry can distinguish between sky and telescope noise, and the 
overscattered tail of the pulse; comparison of the correlated flux density at the peak of the 
pulse, and between pulses, thus provides a measure of the parameter 7 in Figures 2 and 3. 

We are proceeding with additional observational tests that can distinguish among mod- 
els for scattering. As Figures 2,3 suggest, different values of 7 yield different delays between 
the arrival time for the unscattered pulse and the peak of the pulse, as well as different 
behavior in the tail. Observations at short wavelengths yield the arrival time for the un- 
scattered pulse, with corrections for dispersion in the uniform interstellar plasma. Thus, 
multi-frequency observations of pulse broadening, with timing information, can distinguish 
among the different models. 

Shapes of scattered images provide another possible observational test of our theory. 
As Equation (9) suggests, the distribution of deflections of waves, at a given point, is drawn 
from a Levy distribution. Because the Levy distribution is a stable attractor, iteration of the 
process yields a Levy distribution of directions, at any distance from the source. Compared 
with a Gaussian distribution of similar width, a Levy distribution has excesses at both large 
and small angles: it produces large deflections more commonly, and compensates with a 
large population with small deflection. Thus, the scattered image for a Levy flight should 
have a halo at large angles, and a cusp at small angles, compared with a Gaussian image. 
Interestingly, observations of some heavily-scattered extragalactic sources show such core- 
halo excesses (Desai & Fey 2001). Intrinsic structure of the source could possibly explain 
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the observed halo structure; but cannot explain the relatively less-scattered cusp. 



6. 



Discussions and conclusions. 



In the present paper, we have investigated the observational pulse shape. Another object 
directly related to the electron distribution function is the visibility function, ($(x 1 , £)<&*(x 2 , t)}, 
measured at the locations x x and x 2 at the Earth. However, the electron fluctuations probed 
by this function are limited to small scales, |xi — x 2 | < 10 8 cm, which can be comparable 
to the inner scale of the turbulence (Spangler & Gwinn 1990; Kaspi & Stinebring 1992; 
Gupta, Rickett & Coles 1993) where the statistics of the ISM are different from the statistics 
responsible for pulse broadening. The visibility function can be easily constructed in our 
multi-screen model by solving Eqs. (4,6) for AA; = 0, however, we leave its analysis for future 
communication. Here, we estimate what restrictions would be imposed on our theory by an 
inner scale of the order y in ~ 10 8 cm. 

Due to the formulas of Section 3, large-angle refraction is provided by small-scale fluc- 
tuations of electron density. In the case of Levy statistics, a sum of random variables can 
be dominated by a single term, therefore, from formula (17) we can estimate (AC, max ) 2 ~ 
2-ktc/(z N\ 2 B 2 ). From formula (15), the electron-density fluctuations producing such an 
angle deviation should be strong at the scale y min ~ 2n/(BA^ max ) ~ [2A 2 iV ~z / '(tc)] 1 / 2 , or 



where we introduced the dispersion measure DM = (n)Nz , and the signal frequency, v = 
cj A. Quite conveniently, formula (22) includes only the quantities averaged along the line of 
sight, but not the particular parameters I, l , An , and a of the scattering screens. 

As an example, let us consider the parameters for the large-dispersion-measure pulsar 
J1852+0031 (DM=680.0) observed at 1175 MHz in (Bhat et al 2004). Substituting v ~ 
10 9 s _1 , DM ~ 2 • 10 21 cm~ 2 , r ~ 0.5 s, and assuming that (n) ~ 3 • 10~ 2 cm~ 3 , we obtain 
Drain ~ 10 8 cm. Therefore, the condition y min ~ y in is marginally satisfied. Due to the 
scaling r oc \ A DM A , smaller-dispersion-measure pulsars satisfy this condition even better. 
For the cases when this condition is violated, a more appropriate description will be given 
by "truncated" Levy distributions, see e.g., (Nakao 2000). 

In conclusion we would like to mention the observational correlation functions that are 
expressed through the higher-than-second moments of the wave amplitude, $(x, i). These 
observational quantities require a special consideration that at present is not available to 




(22) 



-19- 



us. The examples of such quantities are given by the frequency correlation function that 
is used to characterize pulse broadening for DM < 100 pc cm -3 , see, e.g., (Gwinn et al 
1998), or by the intensity correlation functions that are used to characterize the density 
statistics on scales larger than the Earth size (Shishov et al 2003). The scaling of such 
functions may be easily found in the Gaussian theory when these functions can be related 
to the second-order moment, ((An) 2 ). However, the meaning of these correlation functions 
in the Levy case has yet to be determined. For instance, it seems possible that the deviation 
of the pulse-width scaling from r oc DM 4 for the lower part of the elbow diagram, in 
Fig. (1), is due to the different method of reconstructing r in this region. Indeed, below 
DM ~ lOOpc cm -3 , r is not measured directly, but rather is defined as the inverse frequency 
decorrelation bandwidth (Sutton 1971; Rickett 1977, 1990; Ramachandran, et al 1997; Gwinn 
et al 1998; Bhat et al 2004). 

To summarize, in the present paper we proposed a multi-screen model of the uniform, 
non-Gaussian interstellar medium, and we found its analytic solution. We demonstrate that 
this model explains both the shapes and the scalings of the observational pulse profiles, and it 
is free of the Sutton and Williamson paradoxes that could hardly be resolved in the standard 
Gaussian theory of scintillations. Also, we provided a physical model of density distribution 
having random discontinuities, which produces Levy scintillations and which may serve as a 
good approximation for the analysis of the observational data. We have demonstrated that 
the observational pulse shapes depend on the entire shape of the distribution of electron 
density, rather than simply on the second moment. Pulse profiles may therefore serve as a 
valuable tool for reconstructing this function from observations, and we have developed the 
corresponding method in the present paper. 
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Dispersion Measure DM (cm -3 pc) 



Fig. 1. — Pulse temporal broadening as a function of the dispersion measure, DM = 
J^n(z)dz, which is a measure of the distance to the pulsar. Except as noted, data were 
taken from (Phillips & Clegg 1992). The solid line has slope 4, which contradicts the stan- 
gard Gaussian model of scintillations predicting slope 2 (the Sutton paradox), and which 
agrees with the Levy model with (3 = 2/3. 
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Fig. 2. — Intensities of pulses predicted by the Gaussian model (7 = 2), by the Levy-Cauchy 
model (7 = 1), and by the Levy model with 7 = 2/3, for statistically uniform distribution 
of electron density is the interstellar medium. Note the sharp rise and the pointed apex of 
the signal predicted by the Levy model compared to the signal predicted by the Gaussian 
model. The shapes are rescaled to have similar decaying parts. 
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Fig. 3. — Pulse broadening function for 7 = 2/3 (dashed line), 7=1 (Cauchy: solid line), 
and 7 = 2 (Gaussian: dotted line). The unscattered signal would arrive at time t — 0. Note 
the different relative arrival times of the peaks, as well as the different power-law behavior 
in the tail, for different values of 7. The shapes are rescaled to have same amplitudes and 
half-widths. 
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Fig. 4. — The pulse shape of the pulsar J1848-0123, obtained by Ramachandran, Mitra, 
Deshpande, Connell, & Abies (1997) (Courtesy of Ramachandran Rajagopalan). The data 
are obtained with the Ooty Radio telescope at 327 MHz. Note the sharp rise and the pointed 
apex of the signal. 
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Fig. 5. — Best-fitting pulse-broadening functions for the pulse of pulsar J1848-0123. Solid 
line: 7 = 1 (Cauchy), dotted line: 7 = 2 (Gaussian). Histogram shows the pulse obtained 
by Ramachandran et al. (1997). [Courtesy of Ramachandran Rajagopalan.] Horizontal lines 
at bottom show the zero levels for the two functions; note that the intensity never reaches 
zero for 7 = 1, as a consequence of overscattering. Also note that the Cauchy curve yields 
a steeper rise and a sharper peak, in agreement with the data. 



